Restoring the density-gradient expansion for exchange in solids and surfaces 
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Popular modern generalized gradient approximations (GGA's) are biased toward the description 
of free-atom energies. Restoration of the first-principles gradient expansion for exchange over a 
wide range of density gradients eliminates this bias. We introduce PBEsol, a revised Perdew-Burke- 
Ernzerhof GGA that improves equilibrium properties of densely-packed solids and their surfaces. 
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Ground-state Kohn-Sham density functional theory 
(DFT) [1] has been hugely successful for electronic struc- 
ture calculations of solids and molecules. It relies upon 
good approximations to the exchange-correlation (xc) en- 
ergy as a functional of the electronic spin densities '2'|. 
Over four decades, increasingly sophisticated approxima- 
tions have been developed [3|. The most commonly used 
in solid-state calculations today is the PBE version [4| 
of the generalized gradient approximation (GGA), em- 
ploying both the density and its gradient at each point 
in space. Popular GGAs represent a well-tempered bal- 
ance between computational efficiency, numerical accu- 
racy, and reliability, but PBE also juggles the demands 
of quantum chemistry and solid-state physics [5[. 

While PBE represented a high-point of non-empirical 
functional development 11 years ago, much has since 
been learned about its limitations. PBE reduces the 
chronic overbinding of the local spin density approxi- 
mation (LSDA)[H but, while LSDA often shghtly un- 
derestimates equilibrium lattice constants by about 1%, 
PBE usually overestimates them by about the same 
amount. Other equilibrium properties, such as bulk mod- 
uli, phonon frequencies, magnetism, and ferro-electricity, 
are sensitive to the lattice constant, and so are also over- 
corrected by PBE[g. Surface energies are too low in 
LSDA, but are made lower still by PBE0I. 

However, attempts to construct a better GGA face a 
Procrustean dilemma [Si : Those with an enhanced gra- 
dient dependence [9|, llO| improve atomization and total 
energies, but worsen bond lengths, while more recent 
suggestions of a GGA for solids @, |ll|, [H, [H, El (like 
the very first GGA[153) have a reduced gradient depen- 
dence and typically do improve lattice parameters and/or 
surface energies, but have been criticized for worsen- 
ing total and atomization energies [l3J • More advanced 
functionals have been constructed. For example, meta- 
GGA's, using also the orbital kinetic-energy densities, 
provide greater accuracy over a wider range of systems 
and properties [3|. But current meta-GGA's do not im- 



prove lattice constants as dramatically as surface ener- 
gies, and meta-GGA's are not yet available in all solid- 
state codes. 

In the present work, we explain the origin of this 
dilemma and show that no GGA can do both: Accurate 
atomic exchange energies require violating the gradient 
expansion for slowly- varying densities, which is valid for 
solids and their surfaces. At the GGA level, one must 
choose. A pragmatic approach to lattice properties is 
therefore to use a modified functional especially for solids 
which, unlike previous suggestions, recovers the gradient 
expansion for exchange over a wide range of density gra- 
dients. This by itself yields much better lattice constants 
than PBE (see Ref. [16|), and excellent jellium surface 
exchange energies, because such energies are dominated 
by moderate density gradients. Numerical studies [2| sug- 
gest that PBEsol will also yield highly accurate bulk ex- 
change energies of metals within the pseudopotential ap- 
proximation. These observations suggest an improved 
starting point for more advanced functional construction. 

Restoration of the gradient expanson for exchange re- 
quires a complementary alteration for correlation. We 
fit the jellium surface xc energy, as done by Armiento 
and Mattsson^]. This condition is insufficient to yield 
accurate lattice constants, as experience with the TPSS 
meta-GGA showsfS]- But, as explained below, it is rele- 
vant to lattice constants, because of its bulk-like nature. 
Being accurate for both surface exchange and xc energies, 
PBEsol minimizes reliance on error cancellation between 
X and c, while Ref. [I4I does not. By using the PBE form 
but simply altering two parameters, we retain all other 
exact conditions that make PBE so reliable. 

The GGA form for the exchange energy is simply 
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where n(r) is the electronic density, e™'^(n) is the ex- 
change energy density of a uniform electron gas ( ~ n'^'^), 
s = \\/n\/{2kFn) (with kp = (Sn'^ny/^) is the di- 



mensionless density gradient, and Fx{s) is the enhance- 
ment factor for the given GGAJ^. Eq. (1) is the spin- 
unpolarized form, from which the spin dependence can 
be deduced[4|. Any GGA that recovers the uniform gas 
hmit has 



F,(s) = l + /is2 



(s^O) 



(2) 



The gradient expansion [ll th at is accurate for slowly- 
varying electron gases has[17| /j^ge = 10/81 « .1235. 

To begin, Ref [18| showed that the exchange ener- 
gies of neutral atoms are very well approximated by 
their asymptotic expansion for large Z, i.e., E^ — 
-.2208^5/3 _ iQQz + .... The first term arises from 
LSDA, but the second arises in a GGA from the s^ con- 
tribution to Eq. (2) and requires /i ~ '2^j,ge- Thus any 
GGA that is accurate for the exchange energies of free 
neutral atoms must have fu, « 2/xg_e- PBE does, although 
its value of /z = 0.2195 was found from a different non- 
empirical argument. So does B88, as it was fitted to the 
X energies of noble gas atoms[l9J. Even PW91 does[20|, 
as it reverts only at irrelevantly small s to Hge- 

Thus, to attain accurate exchange energies of atoms 
(vital to dissociation energies in molecules and cohesive 
energies in solids), any GGA must strongly violate the 
gradient expansion for slowly- varying densities [l8|. But 
most of thermochemistry occurs without free atoms, and 
is not much worse in LSDA than in PBE (e.g., [2l|). 
Moreover, for the evaluation of exchange, the densities 
of real solids and their surfaces are often almost slowly- 
varying over space. Restoring the gradient expansion 
should improve their description (but worsen atomization 
energies). The GGA is a limited form, and cannot satisfy 
both conditions. Eq. (2) suggests a necessary condition 
for convergence of the second-order gradient expansion 
for exchange: s ^ 1. Since s ^ 1 for valence electrons 
in densely-packed solids (or s < 2 in core- valence regions 
of alkali atoms) , and since the reduced Laplacian of the 
density is also ^ 1, the gradient expansion is important 
for exchange in solids. We choose ^ge for PBEsol. 

Now, for a GGA correlation functional that recovers 
the uniform gas limit, the gradient expansion is 



Ec[n]^ f d 
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where e"'"^(n) is the correlation energy per particle of the 
uniform gas, /3 is a coefficient, and t = | Vn|/{2 kxp n} is 
the appropriate reduced density gradient for correlation 
(fixed by the Thomas-Fermi screening wavevector fcr f = 
-\/4fcF/7r, not kp-) For slowly- varying high densities [22], 
Pge = 0.0667. Unlike exchange, the second-order term 
in the gradient expansion for correlation cannot be small 
compared to the local term everywhere even for valence 
electrons in solids: /St^ can be large compared to |e""'^| 
(as PGEt'^ = O.ls^/r-s). The gradient expansion can be 
relevant to real systems (especially solids) for exchange, 
but much less so for correlation. 



Also relevant to solids is fxc{q) for the response of the 
uniform gas to a weak potential A cos(q- r). For valence- 
electron densities, the exact fxcil) is almost independent 
of q, up to 2kp [23]. Thus LSDA, which produces a con- 
stant (the value at g = 0), yields an accurate approx- 
imation for q < 2kF- But any GGA with a non-zero 
xc contribution to second order in Vn produces a term 
quadratic in q. Since we are interested in weakly-varying 
valence electron densities in densely-packed solids, we 
wish to retain this excellent feature of LSDA. If 



^i = ^2^/3, 



(4) 



there is complete cancellation between beyond-LSDA 
X and c contributions, restoring LSDA response. 

In PBE, the gradient expansion for correlation is re- 
spected, i.e., P = Pge, and /i « 2^ge satisfies Eq. ([4]). 
This choice agrees well with the PW91 exchange func- 
tional, and with its hole construction [20[ (where sharp 
real-space cutoffs are most appropriate to atoms), and 
yields highly accurate exchange energies of atoms. But 
we have already argued that ^ w 2^ge is harmful for 
many condensed matter applications. Once we choose 
fiGE for exchange, we cannot recover simultaneously the 
GEA for correlation and the linear response of LSDA for 
a uniform density. Exact satisfaction of Eq. ((4]) would 
yield /3 = 0.0375, but a compromise value will satisfy 
another, more relevant constraint for solid-state applica- 
tions. 

For correlation, large neutral jellium clusters are our 



paradigm, for which E^c 
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radius grows, where e|"' is the jellium surface xc energy, 
V the volume of the cluster and A its area. A GGA that 
recovers e|"''^ will be correct in leading- and next-order 
for neutral jellium clusters as A^ — * oo, in a similar way to 
popular exchange GGA's for neutral atoms. Moreover, 
the surface energy is dominated by xc contributions and 
e™''^ is a bulk-like property, arising mainly (103% at r^ = 
2) from a moderately-varying-density region (with s ^ 
1.9) inside the classical turning plane. 

We check that this condition is compatible with the 
restoration of the gradient expansion for exchange. Be- 
cause jellium clusters have a uniform bulk density and 
because most of the surface energy comes from within, 
the gradient expansion should be accurate. We find, at 
bulk density rg = 3, the errors of the surface exchange 
energy are: LSDA 27%, PBE -11% and PBEsol 2.7%. 

We fit e|"''^[12] to determine our correlation functional. 
But the jellium e™''^ is not known exactly. Figure 1 
shows surface energy enhancements relative to LSDA. 
The likely "range of the possible" fo r e™ ''^ extends from 
TPSS meta-GGA [1, or RPA-h [IJ at the low end 
of what is possible (in agreement with the most recent 
Quantum Monte Carlo calculations |25|]) to the RPA-like 
Pitarke-Perdew (PP)[2a] value at the high end. For a 
history of this subject, and another sophisticated esti- 
mate, see Ref. [27|. TPSS may provide the best target 



FIG. 1: Ratio of calculated surface exchange-correlation to 
that of LSDA as a function of rs for various approximations. 



for a GGA. We choose f3 — 0.046 and fi — fioE (within 
the PBE form) for PBEsol, to best fit the TPSS results. 
PBEsol should improve most surface energies over LSDA, 
whereas PBE worsens them. 

Thus we have violated Eq. (|4]) in favor of good sur- 
face energies. But our value for /3 is considerably closer 
to that of the linear response requirement (0.0375) than 
that demanded by complete restoration of the gradient 
expansion (0.0667). 

PBEsol becomes exact for solids under intense com- 
pression, where real solids and their surfaces become 
truly slowly- varying, and exchange dominates over 
corr elation [l3|. In Fig. 2, we plot the enhancement fac- 
tors of PBE and PBEsol. For a spin-unpolarized (C = 0) 
density n — 3/(47rrg), we define F^drs, s) by 



Eg^^[n] = / dVer'(n(r))F.e(rs(r),s(r)). (5) 



The high-density (rs -^ 0) limit is F^{s) of Eq. (1). The 
nonlocality or s-dependence of GGA exchange is dimin- 
ished from PBE to PBEsol, making the latter somewhat 
closer to LSDA. Over the whole range s ^ 1, the PBEsol 
Fx is close to 1 + iices'^- The range < s < 3 is energeti- 
cally important for most properties of most real systems, 



while 



< 



< 



1 and 1 < rs < 10 are the ranges for 



valence-electron regions in many densely-packed solids. 

To test our functional, we employ a test set of 18 solids 
from Ref [28]. These come in four groups: Simple metals 
(Li,Na,K,Al), semiconductors (C,Si,SiC,Ge,GaAs), ionic 
solids (NaF,NaCl,LiCl,LiF,MgO), and transition metals 
(Cu,Rh,Pd,Ag). The set is not claimed to be represen- 
tative, but was chosen for the availability of basis func- 
tions and anharmonic corrections [28j. Our calculations 
use the Gaussian orbital periodic code of Ref 28], with 
basis sets of the same or higher quality. In Table I, wc 
list both mean and mean absolute errors for lattice con- 
stants. The systematic PBE overestimate is close to the 
systematic LSDA underestimate, as shown by the total 



FIG. 2: Enhancement factors of PBE and PBEsol, for spin- 
unpolarized systems, as a function of reduced density gradi- 
ent, for various values of rs. 



class 


LSDA PBE TPSS PBEsol 


mean error 






4 simple metals 


-9.0 


2.9 


5.3 


-0.3 


5 semiconductors 


-1.1 


7.9 


6.2 


3.0 


5 ionic solids 


-8.4 


8.5 


6.8 


2.0 


4 transition metals 


-4.0 


6.4 


2.5 


0.0 


total 


-5.5 


6.6 


5.4 


1.3 


mean absolute error 




4 simple metals 


9.0 


3.4 


5.3 


2.3 


5 semiconductors 


1.3 


7.9 


6.2 


3.0 


5 ionic solids 


8.4 


8.5 


6.8 


2.7 


4 transition metals 


4.0 


6.4 


2.7 


1.9 


total 


5.6 


6.7 


5.4 


2.5 



TABLE I: Errors in equilibrium lattice constants (in Ax 10^^) 
on our data set of 18 solids, relative to experiment with esti- 
mates of the zero-point anharmonic expansion removed [2a]. 



mean absolute errors, and TPSS cures this very little. 
On the other hand, PBEsol greatly reduces this overesti- 
mate, by a factor of almost 4, except for semiconductors, 
where LSDA is unsurpassed. 

PBEsol is not expected to give good atomization ener- 



error 


LSDA PBE TPSS PBEsol 


mean error 
mean abs. error 


3.35 0.54 0.18 1.56 
3.35 0.67 0.26 1.56 



TABLE II; Errors in atomization energies (eV) for the AE6 
set of molecules, using the 6-31H-G(3df,2p) basis set. 



gies. In Table 2, we give the errors on the AEG data set 
of molecules. These 6 molecules (SiH4, S2, SiO, C3H4 
(propyne), C2H2O2 (glyoxal), and C4H8 (cyclobutane) 
were chosen[29| to be representative, i.e., to reproduce 
the errors of much larger data sets. As is clear, and ex- 
pected, PBEsol is much less accurate than PBE, only 
about halving the error of LSDA. This can be related to 
PBEsol's worsened total energies of atoms. 

We have demonstrated the relevance of the second- 
order gradient coefficient for the exchange energy of a 
slowly- varying density to the bulk and surface properties 
of solids. The TPSS meta-GGA [3|, which incorporates 
this coefficient, gets good surface energies but its lattice 
constants are only marginally better than those of PBE 
on which it builds, whereas PBEsol is significantly bet- 
ter. This suggests that an improved meta-GGA needs to 
recover the gradient expansion for exchange over a wider 
range of density distributions n(r) than TPSS does. 

Previous attempts to improve on PBE within the GGA 
form have retained the PBE gradient coefficients^, and B 
for small s, but altered the behavior at large s[6l.l9l.lld|. 



or have zeroed out fj,\lll |12| |. and are thus fundamen- 
tally different from PBEsol. The AM05[il| functional 
performs very similarly to PBEsol for the solids studied 
here, but AMOS follows the proposal of Vitos et al.[ll| to 
fit the conventional exchange energy density of an Airy 
gas. For Ts = 2 to 6, AMOS exchange has errors of 12 
to 60% for el^f, compared to 1.6 to 4.1% for PBEsol. 
Numerical comparisons and details are available [l6|. 

We have identified the simple exchange-correlation 
physics underlying many properties of many solids, and 
shown how it differs from that for atoms. We recommend 
PBEsol for the applications discussed here. Any exist- 
ing code that implements PBE can be instantly mod- 
ified to try PBEsol, by simply replacing the values of 
/i and 13. Modified PBE subroutines are available from 
http://dft.uci.edu. We thank NSF (CHE-03554G5, CHE- 
0457030, and DMR-0501588) and OTKA for support. 
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ERRATUM 

Restoring the density-gradient 

expansion for exchange in sohds and 

surfaces 



The errors of the density-functional lattice constants 
were deviations from corrected experimental values, 
where the corrections removed the finite-temperature 
and zero-point anharmonic expansion (ZPAE) effects. 
The ZPAE corrections were taken from Ref. 28 and 
based on Eq. (18) of Ref. 28, i.e., on Appendix A of 
A.B. Alchagirov et al., Phys. Rev. B 63, 224115 (2001), 
in which e is the energy per atom and vq is the volume 
per atom. But Ref. 28 inadvertently used the volume 
per unit cell, which is twice the volume per atom for 
the semiconductors and ionic crystals. Thus the ZPAE 
corrections for these solids in Ref. 28 and in our article 
were only half what they should have been. Parts of 
Table I (with errors in units of 0.01 A) of our article 
should be corrected as follows: 



class 


LSDA PBE TPSS PBEsol 


mean error 




5 semiconductors 


-0.3 8.6 7.0 


3.7 


5 ionic solids 


-6.9 10.1 8.3 


3.5 


total 


-4.9 7.3 6.0 


1.9 


mean absolute error 




5 semiconductors 


0.9 8.6 7.0 


3.7 


5 ionic solids 


6.9 10.1 8.3 


3.5 


total 


5.1 7.4 6.0 


2.9 



The corrected experimental lattice constants for the 
4 simple metals and 4 transition metals are unaffected. 
No conclusion is changed by this erratum, although the 
improvement of PBEsol relative to LSDA is slightly re- 
duced. Some lattice constants (in A) from the "Expt.- 
ZPAE" column of Table SI of the Supplementary In- 
formation (EPAPS PRLTAO-100-016814) should be cor- 
rected as follows: 3.544 (C), 5.416 (Si), 4.340 (SiC), 5.640 
(Ge), 5.638 (GaAs), 5.566 (NaCl), 4.579 (NaF), 5.074 
(LiCl), 3.964 (LiF), and 4.188 (MgO). J.P.P. thanks Fa- 
bien Tran for pointing out this error. 



